Random restoration projects

Summary plots

# total random sites to create
tot <- nrow(restdat)

id <- stri_rand_strings(tot, 4)
dts <- range(restdat$date)

ext <- bbox(tbpoly)

lon <- rnorm(10 * tot, mean(reststat$lon), sd(reststat$lon))
lat <- rnorm(10 * tot, mean(reststat$lat), sd(reststat$lat))
tmp <- SpatialPoints(cbind(lon, lat), 
                     proj4string = crs(tbpoly)
) %>% 
  .[tbpoly, ] %>% 
  .[sample(1:nrow(.@coords), tot, replace = F), ]

restdat_rnd <- tibble(id) %>% 
  mutate(
    date = sample(seq(dts[1], dts[2]), tot, replace = T),
    top = sample(c('hab', 'wtr'), tot, replace = T)    
  )

reststat_rnd <- tibble(id) %>% 
  mutate(
    lat = tmp$lat, 
    lon = tmp$lon
  )

resgrp <- 'top'
restall_rnd <- left_join(restdat_rnd, reststat_rnd, by = 'id')
names(restall_rnd)[names(restall_rnd) %in% resgrp] <- 'Restoration\ngroup'

# base map
ext <- make_bbox(reststat_rnd$lon, reststat_rnd$lat, f = 0.1)
map <- get_stamenmap(ext, zoom = 10, maptype = "toner-lite")
pbase <- ggmap(map) +
  theme_bw() +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )

# map by restoration type
pbase +
  geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = `Restoration\ngroup`), size = 4, pch = 21)

# map by date
pbase +
  geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = factor(date)), size = 4, pch = 21)

# barplot of date counts
toplo <- restall_rnd %>% 
  group_by(date)
ggplot(restall_rnd, aes(x = factor(date))) + 
  geom_bar() + 
  coord_flip() + 
  theme_bw() + 
  theme(
    axis.title.y = element_blank()
  ) +
  scale_y_discrete(expand = c(0, 0))

Distance to restoration sites

wqmtch_rnd <- get_clo(restdat_rnd, reststat_rnd, wqstat, resgrp = 'top', mtch = mtch)
save(wqmtch_rnd, file = 'data/wqmtch_rnd.RData', compress = 'xz')

head(wqmtch_rnd)

## # A tibble: 6 x 5
##    stat resgrp   rnk    id     dist
##   <int>  <chr> <int> <chr>    <dbl>
## 1    47    hab     1  2L1a 3554.962
## 2    47    hab     2  79vS 4535.917
## 3    47    hab     3  gHCg 6577.988
## 4    47    hab     4  NriW 7685.027
## 5    47    hab     5  osYc 7695.666
## 6    47    hab     6  JFiC 8135.616

Closest

## 
# plots

# combine lat/lon for the plot
toplo <- wqmtch_rnd %>% 
  left_join(wqstat, by = 'stat') %>% 
  left_join(reststat_rnd, by = 'id') %>% 
  rename(
    `Restoration\ngroup` = resgrp,
    `Distance (dd)` = dist
  )
    
# restoration project grouping column
resgrp <- 'top'
restall_rnd <- left_join(restdat_rnd, reststat_rnd, by = 'id')
names(restall_rnd)[names(restall_rnd) %in% resgrp] <- 'Restoration\ngroup'

# extent
ext <- make_bbox(wqstat$lon, wqstat$lat, f = 0.1)
map <- get_stamenmap(ext, zoom = 12, maptype = "toner-lite")

# base map
pbase <- ggmap(map) +
  theme_bw() +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  ) +
  geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = `Restoration\ngroup`), size = 4, pch = 21) +
  geom_point(data = wqstat, aes(x = lon, y = lat), size = 2)

# closest
toplo1 <- filter(toplo, rnk %in% 1)

pbase + 
  geom_segment(data = toplo1, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)

Closest twenty percent

# closest five percent
fvper <- max(toplo$rnk) %>% 
  `*`(0.2) %>% 
  ceiling
toplo2 <- filter(toplo, rnk %in% c(1:fvper))

pbase + 
  geom_segment(data = toplo2, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)

Closest all combinations

# closest all combo
toplo3 <- toplo

pbase + 
  geom_segment(data = toplo3, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)

Summarizing effects of restoration projects

Get weighted average of project type, treatment (before, after) of salinity for all wq station, restoration site combinations.

salchgout_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'sal', yrdf = yrdf, chgout = TRUE)
salchg_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'sal', yrdf = yrdf)
save(salchgout_rnd, file = 'data/salchgout_rnd.RData')
save(salchg_rnd, file = 'data/salchg_rnd.RData')
head(salchgout_rnd)

## # A tibble: 6 x 3
## # Groups:   stat [2]
##    stat     cmb     cval
##   <int>   <chr>    <dbl>
## 1     6 hab_aft 23.90415
## 2     6 hab_bef 23.92485
## 3     6 wtr_aft 23.32146
## 4     6 wtr_bef 23.65890
## 5     7 hab_aft 25.00674
## 6     7 hab_bef 24.19617

head(salchg_rnd)

## # A tibble: 6 x 4
##    stat     hab     wtr     cval
##   <int>  <fctr>  <fctr>    <dbl>
## 1     6 hab_aft wtr_aft 23.61280
## 2     6 hab_aft wtr_bef 23.78152
## 3     6 hab_bef wtr_aft 23.62316
## 4     6 hab_bef wtr_bef 23.79187
## 5     7 hab_aft wtr_aft 24.39071
## 6     7 hab_aft wtr_bef 24.49705

Get conditional probability distributions for the restoration type, treatment effects, salinity as first child node in network.

wqcdt_rnd <- get_cdt(salchg_rnd, 'hab', 'wtr')
head(wqcdt_rnd)

## # A tibble: 4 x 5
##       hab     wtr              data       crv                    prd
##    <fctr>  <fctr>            <list>    <list>                 <list>
## 1 hab_aft wtr_aft <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 2 hab_aft wtr_bef <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 3 hab_bef wtr_aft <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 4 hab_bef wtr_bef <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>

Discretization of salinity conditional probability distributions:

salbrk_rnd <- get_brk(wqcdt_rnd, qts = c(0.33, 0.66), 'hab', 'wtr')
salbrk_rnd

## # A tibble: 8 x 5
##       hab     wtr      qts       brk  clev
##    <fctr>  <fctr>    <dbl>     <dbl> <dbl>
## 1 hab_aft wtr_aft 26.35948 0.4382055     1
## 2 hab_aft wtr_aft 29.98398 0.8755130     2
## 3 hab_aft wtr_bef 26.50050 0.4343718     1
## 4 hab_aft wtr_bef 30.07727 0.8711653     2
## 5 hab_bef wtr_aft 26.29673 0.4453810     1
## 6 hab_bef wtr_aft 29.89586 0.8771097     2
## 7 hab_bef wtr_bef 26.51146 0.4480668     1
## 8 hab_bef wtr_bef 30.02656 0.8747818     2

A plot showing the breaks:

toplo <- dplyr::select(wqcdt_rnd, -data, -crv) %>% 
  unnest
ggplot(toplo, aes(x = cval, y = cumest)) + 
  geom_line() + 
  geom_segment(data = salbrk_rnd, aes(x = qts, y = 0, xend = qts, yend = brk)) +
  geom_segment(data = salbrk_rnd, aes(x = min(toplo$cval), y = brk, xend = qts, yend = brk)) +
  facet_grid(hab ~ wtr) +
  theme_bw()

Get conditional probability distributions for the restoration type, treatment effects, salinity levels, chlorophyll as second child node in network.

# get chlorophyll changes
chlchgout_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'chla', yrdf = yrdf, chgout = TRUE)
chlchg_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'chla', yrdf = yrdf)
save(chlchgout_rnd, file = 'data/chlchgout_rnd.RData')
save(chlchg_rnd, file = 'data/chlchg_rnd.RData')

# merge with salinity, bet salinity levels
salbrk_rnd <- salbrk_rnd %>% 
  group_by(hab, wtr) %>% 
  nest(.key = 'levs')
allchg_rnd <- full_join(chlchg_rnd, salchg_rnd, by = c('hab', 'wtr', 'stat')) %>% 
  rename(
    salev = cval.y, 
    cval = cval.x
  ) %>% 
  group_by(hab, wtr) %>% 
  nest %>% 
  left_join(salbrk_rnd, by = c('hab', 'wtr')) %>% 
  mutate(
    sallev = pmap(list(data, levs), function(data, levs){
      # browser()
      out <- data %>% 
        mutate(
          saval = salev,
          salev = cut(salev, breaks = c(-Inf, levs$qts, Inf), labels = c('lo', 'md', 'hi')),
          salev = as.character(salev)
        )
      
      return(out)
      
    })
  ) %>% 
  dplyr::select(-data, -levs) %>% 
  unnest
salchg_rnd <- dplyr::select(allchg_rnd, stat, hab, wtr, salev, saval)
save(salchg_rnd, file = 'data/salchg_rnd.RData', compress = 'xz')

chlcdt_rnd <- get_cdt(allchg_rnd, 'hab', 'wtr', 'salev')
save(chlcdt_rnd, file = 'data/chlcdt_rnd.RData', compress = 'xz')

chlbrk_rnd <- get_brk(chlcdt_rnd, c(0.33, 0.66), 'hab', 'wtr', 'salev')
chlbrk_rnd %>% 
  print(n = nrow(.))

## # A tibble: 24 x 6
##        hab     wtr salev       qts       brk  clev
##     <fctr>  <fctr> <chr>     <dbl>     <dbl> <dbl>
##  1 hab_aft wtr_aft    lo 12.744188 0.3744332     1
##  2 hab_aft wtr_aft    lo 17.864984 0.8148562     2
##  3 hab_aft wtr_aft    md  7.963425 0.3240829     1
##  4 hab_aft wtr_aft    md  9.933897 0.7877403     2
##  5 hab_aft wtr_aft    hi  4.525562 0.3260305     1
##  6 hab_aft wtr_aft    hi  5.153925 0.7719709     2
##  7 hab_aft wtr_bef    lo 12.527135 0.3840060     1
##  8 hab_aft wtr_bef    lo 17.170391 0.8135389     2
##  9 hab_aft wtr_bef    md  7.534537 0.3079521     1
## 10 hab_aft wtr_bef    md  9.272015 0.7874837     2
## 11 hab_aft wtr_bef    hi  4.289335 0.3932553     1
## 12 hab_aft wtr_bef    hi  4.994263 0.8402122     2
## 13 hab_bef wtr_aft    lo 13.711891 0.4107129     1
## 14 hab_bef wtr_aft    lo 19.820458 0.8408127     2
## 15 hab_bef wtr_aft    md  7.511585 0.3301260     1
## 16 hab_bef wtr_aft    md  9.877967 0.8154849     2
## 17 hab_bef wtr_aft    hi  4.481702 0.3845379     1
## 18 hab_bef wtr_aft    hi  5.199944 0.8133434     2
## 19 hab_bef wtr_bef    lo 13.302918 0.4449256     1
## 20 hab_bef wtr_bef    lo 19.555237 0.8669764     2
## 21 hab_bef wtr_bef    md  7.110791 0.2890397     1
## 22 hab_bef wtr_bef    md  9.272273 0.7448245     2
## 23 hab_bef wtr_bef    hi  4.134515 0.4030903     1
## 24 hab_bef wtr_bef    hi  4.983973 0.8544857     2

Final combinations long format:

chlbar_rnd <- chlbrk_rnd %>% 
  group_by(hab, wtr, salev) %>% 
  nest %>% 
  mutate(
    data = map(data, function(x){
      
      brk <- x$brk
      out <- data.frame(
        lo = brk[1], md = brk[2] - brk[1], hi = 1 - brk[2]
      )
      
      return(out)
      
    })
  ) %>% 
  unnest %>% 
  gather('chllev', 'chlval', lo:hi) %>% 
  mutate(
    salev = factor(salev, levels = c('lo', 'md', 'hi')),
    chllev = factor(chllev, levels = c('lo', 'md', 'hi'))
  )
save(chlbar_rnd, file = 'data/chlbar_rnd.RData', compress = 'xz')

chlbar_rnd %>% 
  print(n = nrow(.))

## # A tibble: 36 x 5
##        hab     wtr  salev chllev    chlval
##     <fctr>  <fctr> <fctr> <fctr>     <dbl>
##  1 hab_aft wtr_aft     lo     lo 0.3744332
##  2 hab_aft wtr_aft     md     lo 0.3240829
##  3 hab_aft wtr_aft     hi     lo 0.3260305
##  4 hab_aft wtr_bef     lo     lo 0.3840060
##  5 hab_aft wtr_bef     md     lo 0.3079521
##  6 hab_aft wtr_bef     hi     lo 0.3932553
##  7 hab_bef wtr_aft     lo     lo 0.4107129
##  8 hab_bef wtr_aft     md     lo 0.3301260
##  9 hab_bef wtr_aft     hi     lo 0.3845379
## 10 hab_bef wtr_bef     lo     lo 0.4449256
## 11 hab_bef wtr_bef     md     lo 0.2890397
## 12 hab_bef wtr_bef     hi     lo 0.4030903
## 13 hab_aft wtr_aft     lo     md 0.4404230
## 14 hab_aft wtr_aft     md     md 0.4636574
## 15 hab_aft wtr_aft     hi     md 0.4459403
## 16 hab_aft wtr_bef     lo     md 0.4295329
## 17 hab_aft wtr_bef     md     md 0.4795316
## 18 hab_aft wtr_bef     hi     md 0.4469569
## 19 hab_bef wtr_aft     lo     md 0.4300998
## 20 hab_bef wtr_aft     md     md 0.4853589
## 21 hab_bef wtr_aft     hi     md 0.4288055
## 22 hab_bef wtr_bef     lo     md 0.4220507
## 23 hab_bef wtr_bef     md     md 0.4557849
## 24 hab_bef wtr_bef     hi     md 0.4513954
## 25 hab_aft wtr_aft     lo     hi 0.1851438
## 26 hab_aft wtr_aft     md     hi 0.2122597
## 27 hab_aft wtr_aft     hi     hi 0.2280291
## 28 hab_aft wtr_bef     lo     hi 0.1864611
## 29 hab_aft wtr_bef     md     hi 0.2125163
## 30 hab_aft wtr_bef     hi     hi 0.1597878
## 31 hab_bef wtr_aft     lo     hi 0.1591873
## 32 hab_bef wtr_aft     md     hi 0.1845151
## 33 hab_bef wtr_aft     hi     hi 0.1866566
## 34 hab_bef wtr_bef     lo     hi 0.1330236
## 35 hab_bef wtr_bef     md     hi 0.2551755
## 36 hab_bef wtr_bef     hi     hi 0.1455143

Discretesize chlorophyll data, all stations:

# discretize all chl data by breaks 
chlbrk_rnd <- chlbrk_rnd %>% 
  group_by(hab, wtr, salev) %>% 
  nest(.key = 'levs')
allchg_rnd <- allchg_rnd %>% 
  group_by(hab, wtr, salev) %>% 
  nest %>% 
  full_join(chlbrk_rnd, by = c('hab', 'wtr', 'salev')) %>% 
  mutate(
    lev = pmap(list(data, levs), function(data, levs){
      
      out <- data %>% 
        mutate(
          lev = cut(cval, breaks = c(-Inf, levs$qts, Inf), labels = c('lo', 'md', 'hi')),
          lev = as.character(lev)
        )

      return(out)
      
    })
  ) %>% 
  dplyr::select(-data, -levs) %>% 
  unnest %>% 
  rename(
    chlev = lev, 
    chval = cval
    )

save(allchg_rnd, file = 'data/allchg_rnd.RData', compress = 'xz')

A bar plot of splits:

ggplot(chlbar_rnd, aes(x = chllev, y = chlval, group = salev, fill = salev)) +
  geom_bar(stat = 'identity', position = 'dodge') +
  facet_grid(hab ~ wtr) +
  theme_bw()

A plot showing the breaks:

toplo <- dplyr::select(chlcdt_rnd, -data, -crv) %>% 
  unnest %>%
  mutate(
    salev = factor(salev, levels = c('lo', 'md', 'hi'))
  )
chlbrk_rnd <- unnest(chlbrk_rnd)
ggplot(toplo, aes(x = cval, y = cumest, group = salev, colour = salev)) + 
  geom_line() + 
  geom_segment(data = chlbrk_rnd, aes(x = qts, y = 0, xend = qts, yend = brk)) +
  geom_segment(data = chlbrk_rnd, aes(x = min(toplo$cval), y = brk, xend = qts, yend = brk)) +
  facet_grid(hab ~ wtr, scales = 'free_x') +
  theme_bw()